Differential Expression Analysis

1 Introduction

Differential Gene Expression (DGE) refers to the process of identifying and quantifying changes in gene expression levels between different biological conditions, such as healthy vs diseased states, treated vs untreated samples, different developmental stages or distinct cell types. Understanding DGE is fundamental in fields like molecular biology, genetics, and biomedical research, as it provides insights into how genes contribute to various biological processes and phenotypes.

Tumor Necrosis Factor (TNF) is a pro-inflammatory cytokine, crucial for the regulation of the immune system’s healthy inflammatory reactions. However, when in overproduction or dysregulation, TNF can contribute to chronic inflammation implicated in various autoimmune and inflammatory diseases, hence the now widespread use of anti-TNF drugs to treat such conditions, also known as TNF inhibitors (Monaco et al. 2014). Despite the success of anti-TNF agents in treating chronic pathological inflammatory reactions, little is known about their impact on the affected tissues at the transcriptome level. As an endeavor to discover the effects of such drugs on the expression levels of different genes, Karagianni and colleagues applied four different anti-TNF drugs on an established mouse model of inflammatory polyarthritis and collected a large number of independent biological replicates from the synovial tissue of healthy, diseased and treated animals (Karagianni et al. 2019). The dataset that pertains to these experiments is used in this workflow, which is part of the study’s computational analyses for identifying and clustering differentially expressed genes. Expression levels were detected and quantified with DNA microarrays, where detected fluorescence indicates the expression of a specific gene against a reference sample.

2 Exploratory Data Analysis

Exploratory Data Analysis (EDA) is a critical process in data science and statistics that involves analyzing and summarizing the main characteristics of a dataset, often using visual methods. The goal of EDA is to uncover patterns, spot anomalies, test hypotheses, and check assumptions with the help of summary statistics and graphical representations (Morgenthaler 2009).

Code 
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(warn = -1)

suppressPackageStartupMessages({
  library(preprocessCore)
  library(umap)
  library(ggplot2)
  library(multcomp)
  library(gplots)
  library(factoextra)
  library(dplyr)
  library(kableExtra)
  library(gprofiler2)
  library(randomForest)
  library(caret)
  library(cowplot)
  library(RColorBrewer)
  library(plotly)
})
Code 
# Read file
data = read.delim(file = "Raw_common18704genes_antiTNF.tsv",
                  header = T, 
                  row.names = 1,
                  sep = "\t")

# Keep gene and sample names
Gene = rownames(data)
Sample = colnames(data)
Code 
boxplot(data, horizontal = T, las = 1, cex.axis = 0.5)
Figure 1: Boxplot of data before quantile normalization
Code 
# Show head of dataframe
kable(head(data)) |>
  kable_styling(bootstrap_options = c("striped")) |>
  scroll_box(width = "100%", height = "100%") |>
  kable_classic()
Exp1Wt_1 Exp1Wt_2 Exp1Wt_3 Exp2Wt_1 Exp2Wt_2 Exp2Wt_3 Exp2Wt_4 Exp3Wt_1 Exp3Wt_2 Exp3Wt_3 Exp1Tg_1 Exp1Tg_2 Exp1Tg_3 Exp2Tg_1 Exp2Tg_2 Exp2Tg_3 Exp2Tg_4 Exp3Tg_1 Exp3Tg_2 Exp3Tg_3 Exp4Tg_1 Exp4Tg_2 Exp4Tg_3 Exp1Rem3_1 Exp1Rem3_2 Exp1Rem3_3 Exp1Rem6_1 Exp1Rem6_2 Exp1Rem6_3 Exp2Rem6_1 Exp2Rem6_2 Exp2Rem6_3 Exp2Rem6_4 Exp3Rem6_1 Exp3Rem6_2 Exp3Rem6_3 Exp2Hum6_1 Exp2Hum6_2 Exp2Hum6_3 Exp2Hum6_4 Exp3Hum6_1 Exp3Hum6_2 Exp3Hum6_3 Exp4Hum6_1 Exp4Hum6_2 Exp4Hum6_3 Exp2Enb6_1 Exp2Enb6_2 Exp2Enb6_3 Exp2Enb6_4 Exp4Enb6_1 Exp4Enb6_2 Exp4Enb6_3 Exp4Enb6_4 Exp4Enb6_5 Exp4Enb6_6 Exp2Cim6_1 Exp2Cim6_2 Exp2Cim6_3 Exp2Cim6_4 Exp4Cim6_1 Exp4Cim6_2 Exp4Cim6_3 Exp4Cim6_4 Exp4Cim6_5 Exp4Cim6_6
A1bg 3.78911 3.90081 3.72526 4.23290 4.14451 4.17460 4.19904 4.37631 4.35622 4.35167 3.78454 4.07887 4.13286 4.05313 4.20157 4.18684 4.27081 4.50461 4.29096 4.36422 3.74793 3.74200 3.79036 3.59538 3.96438 3.79821 3.55839 3.72886 3.73221 4.00460 4.06560 4.10832 4.14768 4.45705 4.40942 4.29107 4.23619 4.11532 4.04981 4.05459 4.38727 4.36555 4.49019 3.55204 3.45422 3.75162 4.11680 4.11782 4.09583 4.05194 3.70226 3.75635 3.82224 3.39433 3.74648 3.48890 4.11928 4.14522 4.10418 4.10931 3.61480 3.37258 3.58425 3.88982 3.28356 4.04937
A1cf 4.21654 4.17014 4.05211 4.60213 4.62320 4.65975 4.57693 4.15842 4.21211 4.18009 3.98771 4.11709 4.10684 4.68407 4.63325 4.73126 4.67542 4.09641 3.97244 4.10268 3.30531 3.29409 3.06133 4.29460 4.22492 4.08466 4.08205 4.13011 4.09484 4.62164 4.66620 4.61722 4.55999 4.38385 4.06813 4.13370 4.53864 4.75150 4.64363 4.68842 4.03449 4.62612 4.24503 3.39752 3.35704 3.31550 4.61107 4.71722 4.65451 4.68732 3.34686 3.26266 3.04503 3.16430 3.13805 3.10227 4.65152 4.67835 4.67399 4.61035 2.90180 3.27259 3.02554 3.12053 3.08666 3.24268
A2ld1 5.32714 5.62827 5.66799 8.32861 8.27225 8.27447 8.42341 8.22626 8.39943 8.57076 5.23156 5.29504 5.28384 7.90894 7.85233 8.14223 7.79835 7.90913 8.08984 7.95827 6.55196 6.51255 6.43711 6.01249 5.46742 5.54192 5.56968 5.70214 5.46825 8.26541 8.30650 8.22540 8.09633 8.30025 7.86106 7.86827 8.19100 8.37507 8.20061 8.27231 8.39751 8.30546 8.49839 6.69627 6.52055 6.64086 8.00233 8.18752 8.25757 8.03628 6.43991 6.62674 6.80411 6.40195 6.77547 6.49056 8.09591 8.09534 8.08303 8.07122 6.51077 6.56018 6.40339 6.71310 6.35072 6.61334
A2m 5.28545 6.49319 5.19705 6.34696 5.68818 7.02162 6.30722 7.24319 6.80225 6.69186 4.71382 5.03411 5.41249 6.57654 6.10905 6.11674 5.67339 6.93047 6.70562 6.65830 5.29271 5.47738 6.79552 6.19873 5.40119 6.48264 5.78780 5.32115 5.44562 6.88086 5.34267 6.33201 6.39798 6.45636 6.43533 6.19721 6.10152 6.31918 5.95871 6.76260 6.76862 6.14994 5.71077 5.23411 5.38092 5.36692 5.98963 5.63186 5.84172 5.34126 5.51386 5.58334 5.51965 5.63152 5.64672 6.46401 5.42400 6.97807 6.39452 5.24073 4.64791 4.64806 6.10691 5.04173 4.90970 4.32605
A3galt2 4.59905 4.59698 4.52439 4.85959 5.07224 5.01061 5.07838 5.62244 5.79039 5.66514 4.73714 4.90568 5.18567 4.95670 4.65065 4.95141 4.96971 5.23615 5.38050 5.49452 4.72958 4.78397 4.66467 4.41278 4.82166 4.49756 4.54898 4.53977 4.85098 4.89872 4.96186 5.01156 5.02533 5.50526 5.51928 5.52995 5.12319 5.15686 5.02892 4.91145 5.58904 5.62234 5.64090 4.49644 4.98682 4.62059 4.94060 5.10253 4.92040 4.80770 4.63297 4.79730 4.73491 4.63326 4.78748 4.69899 4.96660 5.02830 4.96265 5.03972 4.70114 4.64965 4.50466 4.70696 4.55098 4.64297
A4galt 7.73303 7.65939 7.93984 8.55175 8.48415 8.47890 8.22892 8.72356 8.71276 8.58030 7.22117 7.00438 6.80762 8.65125 8.49504 8.61091 8.60534 8.90632 8.84868 8.88411 8.87908 8.68815 8.59815 7.98948 7.60469 7.47936 7.60934 7.44681 7.36165 8.70350 8.61929 8.80555 8.47461 8.76373 9.06743 9.36881 8.55062 8.73755 8.65360 8.69017 9.22244 9.12786 9.25526 8.76674 8.93546 8.86313 8.66064 8.60910 8.63829 8.86686 8.73614 8.96106 8.77789 8.77312 8.97563 8.70684 8.83717 8.88996 8.98343 8.75343 8.80162 8.76917 8.70345 8.80999 8.86269 8.62107

3 Check for missing values

If there are any missing values, we will need to decide how to handle them, probably by removing the respective genes.

Code 
# Check genes_data for missing values
colSums(is.na(data))
  Exp1Wt_1   Exp1Wt_2   Exp1Wt_3   Exp2Wt_1   Exp2Wt_2   Exp2Wt_3   Exp2Wt_4 
         0          0          0          0          0          0          0 
  Exp3Wt_1   Exp3Wt_2   Exp3Wt_3   Exp1Tg_1   Exp1Tg_2   Exp1Tg_3   Exp2Tg_1 
         0          0          0          0          0          0          0 
  Exp2Tg_2   Exp2Tg_3   Exp2Tg_4   Exp3Tg_1   Exp3Tg_2   Exp3Tg_3   Exp4Tg_1 
         0          0          0          0          0          0          0 
  Exp4Tg_2   Exp4Tg_3 Exp1Rem3_1 Exp1Rem3_2 Exp1Rem3_3 Exp1Rem6_1 Exp1Rem6_2 
         0          0          0          0          0          0          0 
Exp1Rem6_3 Exp2Rem6_1 Exp2Rem6_2 Exp2Rem6_3 Exp2Rem6_4 Exp3Rem6_1 Exp3Rem6_2 
         0          0          0          0          0          0          0 
Exp3Rem6_3 Exp2Hum6_1 Exp2Hum6_2 Exp2Hum6_3 Exp2Hum6_4 Exp3Hum6_1 Exp3Hum6_2 
         0          0          0          0          0          0          0 
Exp3Hum6_3 Exp4Hum6_1 Exp4Hum6_2 Exp4Hum6_3 Exp2Enb6_1 Exp2Enb6_2 Exp2Enb6_3 
         0          0          0          0          0          0          0 
Exp2Enb6_4 Exp4Enb6_1 Exp4Enb6_2 Exp4Enb6_3 Exp4Enb6_4 Exp4Enb6_5 Exp4Enb6_6 
         0          0          0          0          0          0          0 
Exp2Cim6_1 Exp2Cim6_2 Exp2Cim6_3 Exp2Cim6_4 Exp4Cim6_1 Exp4Cim6_2 Exp4Cim6_3 
         0          0          0          0          0          0          0 
Exp4Cim6_4 Exp4Cim6_5 Exp4Cim6_6 
         0          0          0

4 Data Distribution

Data distribution refers to how data values are spread or dispersed across a range of possible values. Understanding the distribution of the data is a key part of Exploratory Data Analysis (EDA), as it gives insights into the central tendency, variability, and shape of the dataset. It can help in detecting patterns, outliers, and important characteristics of the data, and also informs which statistical techniques to use. Thus, a detailed assessment of data distribution is a foundational step in unraveling the complexities of any dataset during the EDA process.

Central Tendency

  • Mean (Average): The sum of all values divided by the number of values. It provides a sense of the center of the data.

  • Median: The middle value when data is ordered. It’s less affected by extreme values (outliers) and provides a robust measure of the center.

  • Mode: The most frequent value in the data. For categorical data, this is often the only measure of central tendency.

Spread (Dispersion)

  • Range: The difference between the maximum and minimum values.

  • Variance: The average of the squared differences from the mean. It quantifies how much the data points deviate from the mean.

  • Standard Deviation: The square root of the variance. It provides a more interpretable measure of spread in the same units as the data.

  • Interquartile Range (IQR): The range within the middle 50% of the data (between the first and third quartiles). It’s a robust measure of spread that isn’t affected by outliers.

Code 
# Adjust the layout and margins
par(mfrow = c(8, 9), mar = c(1, 1, 2.5, 1))

for (col in colnames(data)) {
  plot(density(data[[col]]), 
       main = col,
       xlab = col, col = "#009AEF", lwd = 2)
}
Figure 2: Data distribution before quantile normalization
Code 
data %>% 
  select_if(is.numeric) %>%
  apply(2, function(x) round(summary(x), 3)) %>% 
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "bordered")) %>% 
  kable_classic() %>%
  scroll_box(width = "100%", height = "100%")
Exp1Wt_1 Exp1Wt_2 Exp1Wt_3 Exp2Wt_1 Exp2Wt_2 Exp2Wt_3 Exp2Wt_4 Exp3Wt_1 Exp3Wt_2 Exp3Wt_3 Exp1Tg_1 Exp1Tg_2 Exp1Tg_3 Exp2Tg_1 Exp2Tg_2 Exp2Tg_3 Exp2Tg_4 Exp3Tg_1 Exp3Tg_2 Exp3Tg_3 Exp4Tg_1 Exp4Tg_2 Exp4Tg_3 Exp1Rem3_1 Exp1Rem3_2 Exp1Rem3_3 Exp1Rem6_1 Exp1Rem6_2 Exp1Rem6_3 Exp2Rem6_1 Exp2Rem6_2 Exp2Rem6_3 Exp2Rem6_4 Exp3Rem6_1 Exp3Rem6_2 Exp3Rem6_3 Exp2Hum6_1 Exp2Hum6_2 Exp2Hum6_3 Exp2Hum6_4 Exp3Hum6_1 Exp3Hum6_2 Exp3Hum6_3 Exp4Hum6_1 Exp4Hum6_2 Exp4Hum6_3 Exp2Enb6_1 Exp2Enb6_2 Exp2Enb6_3 Exp2Enb6_4 Exp4Enb6_1 Exp4Enb6_2 Exp4Enb6_3 Exp4Enb6_4 Exp4Enb6_5 Exp4Enb6_6 Exp2Cim6_1 Exp2Cim6_2 Exp2Cim6_3 Exp2Cim6_4 Exp4Cim6_1 Exp4Cim6_2 Exp4Cim6_3 Exp4Cim6_4 Exp4Cim6_5 Exp4Cim6_6
Min. 3.376 3.144 3.338 3.726 3.615 3.586 3.555 1.029 0.978 0.806 3.349 3.261 3.282 3.607 3.728 3.754 3.605 0.751 1.088 1.023 1.177 1.101 1.058 3.407 3.372 3.283 3.390 3.329 3.305 3.682 3.624 3.589 3.662 0.978 1.015 1.059 3.626 3.578 3.658 3.689 1.100 0.953 1.021 1.084 1.144 1.114 3.733 3.609 3.716 3.787 1.210 1.096 1.147 1.376 1.283 1.318 3.632 3.676 3.691 3.685 1.233 1.147 1.014 1.377 1.201 1.128
1st Qu. 4.550 4.524 4.441 5.039 4.981 4.993 5.038 5.169 5.199 5.276 4.596 4.719 4.808 4.973 4.938 4.987 5.027 5.155 5.186 5.172 4.358 4.395 4.352 4.565 4.675 4.493 4.525 4.500 4.705 4.952 4.956 4.982 5.004 5.269 5.199 5.177 4.976 4.975 4.952 4.968 5.170 5.201 5.208 4.376 4.365 4.368 5.027 5.018 4.979 4.960 4.366 4.379 4.305 4.365 4.314 4.378 4.973 4.994 4.962 4.976 4.323 4.316 4.282 4.319 4.344 4.301
Median 6.030 6.070 6.107 6.790 6.841 6.833 6.752 6.727 6.706 6.665 6.042 5.997 6.029 6.831 6.840 6.852 6.781 6.734 6.730 6.745 5.806 5.761 5.783 6.126 6.066 6.057 6.040 6.026 6.016 6.813 6.812 6.810 6.849 6.691 6.757 6.714 6.837 6.804 6.827 6.839 6.744 6.726 6.708 5.806 5.805 5.789 6.811 6.808 6.825 6.864 5.765 5.780 5.798 5.773 5.778 5.792 6.834 6.798 6.847 6.815 5.786 5.805 5.799 5.788 5.806 5.762
Mean 6.260 6.285 6.286 6.960 6.952 6.944 6.934 6.623 6.619 6.616 6.270 6.221 6.233 6.966 6.952 6.981 6.943 6.616 6.624 6.626 5.743 5.726 5.730 6.341 6.296 6.266 6.265 6.240 6.245 6.948 6.952 6.952 6.952 6.624 6.637 6.623 6.958 6.958 6.953 6.969 6.624 6.635 6.630 5.756 5.737 5.732 6.957 6.960 6.950 6.972 5.730 5.750 5.738 5.736 5.732 5.744 6.953 6.942 6.957 6.952 5.738 5.740 5.731 5.740 5.747 5.728
3rd Qu. 7.679 7.745 7.797 8.523 8.536 8.526 8.487 8.123 8.092 8.028 7.649 7.452 7.389 8.601 8.594 8.599 8.522 8.147 8.126 8.138 7.113 7.046 7.091 7.826 7.637 7.727 7.704 7.679 7.513 8.556 8.576 8.538 8.524 8.056 8.135 8.133 8.560 8.570 8.568 8.581 8.127 8.131 8.109 7.134 7.091 7.101 8.525 8.532 8.542 8.612 7.068 7.099 7.134 7.094 7.120 7.099 8.555 8.523 8.573 8.558 7.133 7.147 7.140 7.137 7.126 7.129
Max. 13.482 13.700 13.650 13.684 13.776 13.726 13.653 12.926 12.912 12.887 13.204 13.218 13.116 13.725 13.781 13.700 13.700 13.081 12.927 13.040 12.801 12.838 13.148 13.711 13.501 13.562 13.473 13.644 13.453 13.711 13.700 13.694 13.726 12.937 13.001 12.999 13.694 13.699 13.718 13.704 13.073 12.970 12.947 12.823 12.837 12.995 13.695 13.696 13.734 13.767 12.731 12.831 12.784 12.752 12.816 12.797 13.712 13.791 13.753 13.733 12.815 12.854 12.995 12.851 12.833 12.759
Code 
calculate_metrics = function(data.frame) {
  max <- apply(data.frame, 2, max)
  min <- apply(data.frame, 2, min)
  mean <- (max + min) / 2

  dt_matrix = data.frame(name = colnames(data.frame),
                         min = as.numeric(as.character(min)),
                         max = as.numeric(as.character(max)),
                         mean = as.numeric(as.character(mean)))
  return(dt_matrix)
}

# Calculate metrics for each condition
c_metrics = calculate_metrics(data)
c_metrics
         name     min     max     mean
1    Exp1Wt_1 3.37556 13.4822 8.428880
2    Exp1Wt_2 3.14427 13.7005 8.422385
3    Exp1Wt_3 3.33822 13.6498 8.494010
4    Exp2Wt_1 3.72646 13.6836 8.705030
5    Exp2Wt_2 3.61474 13.7756 8.695170
6    Exp2Wt_3 3.58591 13.7262 8.656055
7    Exp2Wt_4 3.55539 13.6532 8.604295
8    Exp3Wt_1 1.02914 12.9258 6.977470
9    Exp3Wt_2 0.97795 12.9121 6.945025
10   Exp3Wt_3 0.80624 12.8875 6.846870
11   Exp1Tg_1 3.34917 13.2042 8.276685
12   Exp1Tg_2 3.26118 13.2179 8.239540
13   Exp1Tg_3 3.28197 13.1161 8.199035
14   Exp2Tg_1 3.60711 13.7250 8.666055
15   Exp2Tg_2 3.72809 13.7809 8.754495
16   Exp2Tg_3 3.75376 13.6998 8.726780
17   Exp2Tg_4 3.60466 13.6999 8.652280
18   Exp3Tg_1 0.75115 13.0814 6.916275
19   Exp3Tg_2 1.08764 12.9265 7.007070
20   Exp3Tg_3 1.02303 13.0399 7.031465
21   Exp4Tg_1 1.17688 12.8012 6.989040
22   Exp4Tg_2 1.10078 12.8384 6.969590
23   Exp4Tg_3 1.05794 13.1482 7.103070
24 Exp1Rem3_1 3.40727 13.7114 8.559335
25 Exp1Rem3_2 3.37224 13.5005 8.436370
26 Exp1Rem3_3 3.28277 13.5621 8.422435
27 Exp1Rem6_1 3.38981 13.4728 8.431305
28 Exp1Rem6_2 3.32883 13.6436 8.486215
29 Exp1Rem6_3 3.30500 13.4533 8.379150
30 Exp2Rem6_1 3.68185 13.7110 8.696425
31 Exp2Rem6_2 3.62371 13.7001 8.661905
32 Exp2Rem6_3 3.58858 13.6942 8.641390
33 Exp2Rem6_4 3.66171 13.7255 8.693605
34 Exp3Rem6_1 0.97805 12.9367 6.957375
35 Exp3Rem6_2 1.01526 13.0005 7.007880
36 Exp3Rem6_3 1.05879 12.9986 7.028695
37 Exp2Hum6_1 3.62648 13.6944 8.660440
38 Exp2Hum6_2 3.57819 13.6989 8.638545
39 Exp2Hum6_3 3.65832 13.7184 8.688360
40 Exp2Hum6_4 3.68866 13.7042 8.696430
41 Exp3Hum6_1 1.09956 13.0734 7.086480
42 Exp3Hum6_2 0.95331 12.9700 6.961655
43 Exp3Hum6_3 1.02125 12.9465 6.983875
44 Exp4Hum6_1 1.08401 12.8233 6.953655
45 Exp4Hum6_2 1.14380 12.8366 6.990200
46 Exp4Hum6_3 1.11394 12.9948 7.054370
47 Exp2Enb6_1 3.73300 13.6948 8.713900
48 Exp2Enb6_2 3.60938 13.6959 8.652640
49 Exp2Enb6_3 3.71600 13.7344 8.725200
50 Exp2Enb6_4 3.78704 13.7671 8.777070
51 Exp4Enb6_1 1.21022 12.7314 6.970810
52 Exp4Enb6_2 1.09627 12.8306 6.963435
53 Exp4Enb6_3 1.14697 12.7839 6.965435
54 Exp4Enb6_4 1.37550 12.7517 7.063600
55 Exp4Enb6_5 1.28300 12.8161 7.049550
56 Exp4Enb6_6 1.31791 12.7969 7.057405
57 Exp2Cim6_1 3.63227 13.7120 8.672135
58 Exp2Cim6_2 3.67638 13.7910 8.733690
59 Exp2Cim6_3 3.69146 13.7532 8.722330
60 Exp2Cim6_4 3.68540 13.7332 8.709300
61 Exp4Cim6_1 1.23322 12.8151 7.024160
62 Exp4Cim6_2 1.14684 12.8535 7.000170
63 Exp4Cim6_3 1.01430 12.9947 7.004500
64 Exp4Cim6_4 1.37744 12.8509 7.114170
65 Exp4Cim6_5 1.20130 12.8331 7.017200
66 Exp4Cim6_6 1.12756 12.7590 6.943280

5 Data Normalization

The primary goal of quantile normalization is to align the statistical distributions of all samples so that they are the same. This technique assumes that all samples should have the same overall distribution of values (e.g., gene expression levels), meaning that any differences between samples that aren’t due to biological variation are artifacts that need to be corrected.

Steps of Quantile Normalization

  • Sort the Data:

    For each sample (e.g., each column in a matrix of gene expression data), the data is sorted / ranked in ascending order.

  • Average the Sorted Values:

    Across all samples, for each rank, the values are averaged. This creates a “reference” distribution, where each rank corresponds to the average value of that rank across all samples.

  • Replace Original Values:

    Each value in the original dataset is replaced by the corresponding value from the reference distribution.

  • Output the Normalized Data:

    The result is a new dataset where each sample has the same distribution of values, which helps remove any technical variations between samples.

Code 
# Convert dataframe to matrix
data = as.matrix(data)

# Normalize data
data = normalize.quantiles(data , copy=TRUE)

# Convert matrix to dataframe
data = data.frame(data)

# Add column names to dataframe
colnames(data) = Sample

# Add row names to dataframe
rownames(data) = Gene
Code 
boxplot(data, horizontal = T, las = 1, cex.axis = 0.5)
Figure 3: Boxplot of data after quantile normalization
Code 
# Adjust the layout and margins 
par(mfrow = c(8, 9), mar = c(1, 1, 2.5, 1))

for (col in colnames(data)) {
  plot(density(data[[col]]), 
       main = col,
       xlab = col, col = "#009AEF", lwd = 2)
}
Figure 4: Data distribution after quantile normalization

6 Dimension Reduction Algorithms

6.1 Uniform Manifold Approximation and Projection (UMAP)

UMAP (Uniform Manifold Approximation and Projection) is a dimensionality reduction algorithm that is widely used in data science and machine learning. It is designed to help visualize high-dimensional data in a lower-dimensional space, typically 2D or 3D, while preserving the global and local structure of the data as much as possible.

Key concepts of UMAP

  1. Manifold Learning: UMAP is based on the idea that high-dimensional data often lies on a lower-dimensional manifold, a geometry concept. UMAP aims to learn this manifold and then project the data onto a lower-dimensional space that best represents the structure of the data.

  2. Topology: UMAP uses topological techniques to understand the shape of the data. It constructs a graph that represents the relationships between data points based on their proximity. This graph is then used to project the data into a lower-dimensional space.

  3. Optimization: The final embedding is obtained by optimizing a cost function that balances the preservation of local structures (how points are grouped) and global structures (the overall layout).

Applications

  1. Visualization: UMAP is often used to visualize complex datasets, such as in biology (e.g., gene expression data), where it helps to identify patterns and clusters.

  2. Preprocessing: It can be used as a preprocessing step before other machine learning algorithms, such as clustering or classification, to reduce dimensionality while maintaining the integrity of the data structure.

  3. Data Exploration: UMAP is useful for exploring and understanding large datasets, making it easier to identify trends, outliers, and relationships within the data.

Comparison to other techniques

  • t-SNE: UMAP is similar to t-SNE, another popular dimensionality reduction technique. However, UMAP is generally faster and better at preserving both global and local structures, making it more effective for large datasets.

  • PCA (Principal Component Analysis): PCA is another common method for dimensionality reduction but focuses on preserving linear relationships, which may not capture the complexity of data as effectively as UMAP, especially in cases of non-linear structures.

Code 
# Keep only WT and TG samples
wt_tg_df = data[, 1:23]

# After dataframe transposition columns must represent genes
wt_tg_df = t(wt_tg_df)
Code 
# UMAP dimension reduction for wt and tg samples
wt_tg_df.umap <- umap(wt_tg_df, n_components=2, random_state=15)

# Keep the numeric dimensions
wt_tg_df.umap <- wt_tg_df.umap[["layout"]]

# Create vector with groups
group = c(rep("A_Wt", 10), rep("B_Tg", 13))

# Create final dataframe with dimensions and group for plotting
wt_tg_df.umap <- cbind(wt_tg_df.umap, group)
wt_tg_df.umap <- data.frame(wt_tg_df.umap)

# Plot UMAP results
ggplotly(
  ggplot(wt_tg_df.umap, aes(x = V1, y = V2, color = group)) +
    geom_point() +
    labs(
      x = "UMAP1",
      y = "UMAP2",
      title = "UMAP plot",
      subtitle = "A UMAP Visualization of WT and TG samples") +
    theme(
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank()
    )
)
Figure 5: UMAP plot

6.2 Principal Component Analysis (PCA)

Principal Component Analysis (PCA) is a statistical technique used to simplify complex datasets by reducing their dimensionality while preserving as much variance as possible. The main idea behind PCA is to transform the original data into a new coordinate system where the greatest variance in the data is captured in the first few dimensions (called principal components). Here’s how PCA works step-by-step:

Covariance Matrix Computation

  • Calculate the covariance matrix of the data. The covariance matrix is a square matrix that shows the covariance between different features. If there are \(n\) features, the covariance matrix will be of size \(n×n\). Covariance gives an idea of how much two features vary together.

Eigenvalue and Eigenvector Calculation

  • Compute the eigenvalues and eigenvectors of the covariance matrix.
  • Eigenvectors represent the directions of the new feature space, and eigenvalues tell us how much variance is along each of these directions.
  • The eigenvector with the highest eigenvalue is the first principal component, which captures the most variance in the data. The second principal component is orthogonal to the first and captures the second most variance, and so on.

Sort Eigenvectors

  • Sort the eigenvectors by their corresponding eigenvalues in descending order. This ordering helps in deciding which principal components to keep. The components corresponding to the largest eigenvalues capture the most significant structure in the data.

Dimensionality Reduction

  • Choose the top \(k\) eigenvectors (where \(k\) is the number of dimensions you want to reduce your data to) and project the original data onto these eigenvectors. This transformation gives you the data represented in terms of the principal components, reducing the dimensionality from \(n\) to \(k\).
  • The resulting \(k\)-dimensional space retains most of the variability in the original \(n\)-dimensional space.

Projection of data

  • Finally, the original data is transformed into the new feature space defined by the selected principal components. This projection is done by multiplying the original data matrix by the matrix of eigenvectors.

Applications

  • Data compression: Reducing the number of features while retaining most of the original information.

  • Noise reduction: By keeping the most significant components, noise in the data (which might be in less significant components) can be reduced.

  • Visualization: Reducing data to 2 or 3 dimensions makes it easier to visualize, especially in high-dimensional datasets.

  • Preprocessing: Used as a preprocessing step before applying machine learning algorithms to improve performance and reduce computational cost.

Code 
# PCA dimension reduction
wt_tg_df.pca <- prcomp(wt_tg_df, scale. = FALSE)
summary(wt_tg_df.pca)
Importance of components:
                           PC1     PC2     PC3      PC4     PC5    PC6     PC7
Standard deviation     64.2329 34.1273 30.6454 20.32780 13.6953 8.9139 7.93867
Proportion of Variance  0.5608  0.1583  0.1277  0.05617  0.0255 0.0108 0.00857
Cumulative Proportion   0.5608  0.7191  0.8468  0.90298  0.9285 0.9393 0.94785
                           PC8     PC9    PC10    PC11    PC12    PC13    PC14
Standard deviation     6.80061 6.34603 6.13851 5.80207 5.57370 5.31339 5.10862
Proportion of Variance 0.00629 0.00547 0.00512 0.00458 0.00422 0.00384 0.00355
Cumulative Proportion  0.95413 0.95961 0.96473 0.96930 0.97353 0.97737 0.98091
                          PC15    PC16    PC17    PC18    PC19    PC20    PC21
Standard deviation     4.95260 4.73197 4.56130 4.14881 3.95569 3.74775 3.64909
Proportion of Variance 0.00333 0.00304 0.00283 0.00234 0.00213 0.00191 0.00181
Cumulative Proportion  0.98425 0.98729 0.99012 0.99246 0.99459 0.99649 0.99830
                         PC22      PC23
Standard deviation     3.5312 8.294e-14
Proportion of Variance 0.0017 0.000e+00
Cumulative Proportion  1.0000 1.000e+00
Code 
plot_grid(fviz_pca_ind(wt_tg_df.pca, repel = TRUE, # Avoid text overlapping
                  habillage = group,
                  label = "none",
                  axes = c(1, 2), # choose PCs to plot
                  addEllipses = TRUE,
                  ellipse.level = 0.95,
                  title = "Biplot: PC1 vs PC2") +
                  scale_color_manual(values = c('#33cc00','#009AEF95')) +
                  scale_fill_manual(values = c('#33cc00','#009AEF95')),
          fviz_pca_ind(wt_tg_df.pca, repel = TRUE, # Avoid text overlapping
                  habillage = group,
                  label = "none",
                  axes = c(1, 3), # choose PCs to plot
                  addEllipses = TRUE,
                  ellipse.level = 0.95,
                  title = "Biplot: PC1 vs PC3") + 
                  scale_color_manual(values = c('#33cc00','#009AEF95')) +
                  scale_fill_manual(values = c('#33cc00','#009AEF95')),
          fviz_pca_ind(wt_tg_df.pca, repel = TRUE, # Avoid text overlapping
                  habillage = group,
                  label = "none",
                  axes = c(2, 3), # choose PCs to plot
                  addEllipses = TRUE,
                  ellipse.level = 0.95,
                  title = "Biplot: PC2 vs PC3") +
                  scale_color_manual(values = c('#33cc00','#009AEF95')) +
                  scale_fill_manual(values = c('#33cc00','#009AEF95')),
          
          # Visualize eigenvalues/variances
          fviz_screeplot(wt_tg_df.pca, 
                    addlabels = TRUE,
                    title = "Principal Components Contribution",
                    ylim = c(0, 65), 
                    barcolor = "#009AEF95",  
                    barfill = "#009AEF95"),
          
          # Contributions of features to PC1
          fviz_contrib(wt_tg_df.pca, 
                  choice = "var", 
                  axes = 1, 
                  top = 14, 
                  color = "#009AEF95", 
                  fill = "#009AEF95"),
          
          # Contributions of features to PC2
          fviz_contrib(wt_tg_df.pca, 
                  choice = "var", 
                  axes = 2, 
                  top = 14, 
                  color = "#009AEF95", 
                  fill = "#009AEF95"),
          labels = c("A", "B", "C", "D", "E", "F")
)
Figure 6: Plot PCA results
Code 
wt_tg_df.pca <- data.frame("PC1" = wt_tg_df.pca$x[,1], 
                           "PC2" = wt_tg_df.pca$x[,2], 
                           "group" = group)

# Plot PCA results

# insert dataframe [1] , variables [2]-[3] and color groyp [4]
ggplotly(
  ggplot( wt_tg_df.pca, aes(x= PC1 , y= PC2 , color= group ))+
    
    # try "geom_point" or "geom_line"
    geom_point()+
    
    # try "ggtitle" or "ggname"
    ggtitle("Two First Components of PCA") +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks = element_blank())
)
Figure 7: Plot two first components of PCA

7 Statistical Analysis

7.1 Group treatments in dataframe

ANOVA, or Analysis of Variance, is a statistical method used to compare means among three or more groups to determine if there are any statistically significant differences between them. It extends the t-test, which is typically used to compare the means of two groups, to situations where more than two groups are involved. ANOVA allows for the simultaneous assessment of variations within and between groups, enabling the identification of genes with expression patterns that are significantly different across experimental conditions. It is particularly useful in experiments where multiple groups or conditions are being tested simultaneously.

Hypotheses in ANOVA

ANOVA tests the following hypotheses:

  • Null Hypothesis (H₀): All group means are equal. (No significant difference between groups)

  • Alternative Hypothesis (H₁): At least one group mean is different from the others. (There is a significant difference between groups)

Partitioning Variance

ANOVA works by partitioning the total variance observed in the data into two components:

  1. Between-Group Variance: The variation due to the interaction between the different groups. This measures how much the group means differ from the overall mean.

  2. Within-Group Variance (Error Variance): The variation within each group. This measures how much individual observations differ from their group mean.

Mathematically, the total variance is expressed as:

$$ \text{Total Variance} = \text{Between-Group Variance} + \text{Within-Group Variance} $$

F-Ratio Calculation

ANOVA calculates an F-ratio by comparing the between-group variance to the within-group variance:

\[ F = \frac{\text{Between-Group Variance}}{\text{Within-Group Variance}} \]

  • A large F-ratio indicates that the between-group variance is large relative to the within-group variance, suggesting that the group means are significantly different.

  • A small F-ratio suggests that any observed differences in means are more likely due to random chance.

ANOVA table

The results of ANOVA are typically summarized in an ANOVA table, which includes:

  • Sum of Squares (SS): Measures of variation for both between-group and within-group components.

  • Degrees of Freedom (df): The number of values that are free to vary for each component.

  • Mean Squares (MS): Calculated by dividing the sum of squares by the corresponding degrees of freedom.

  • F-Statistic (F): The ratio of mean squares between groups to mean squares within groups.

Assumptions of ANOVA

For the results of ANOVA to be valid, several assumptions need to be met:

  1. Independence: The observations within each group and between groups should be independent of each other.

  2. Normality: The data in each group should be approximately normally distributed.

  3. Homogeneity of Variances: The variances among the groups should be approximately equal.

Code 
# -------------------Apply ANOVA on all genes----------------------

# Create Matrix by Excluding rownames and colnames
matrixdata = as.matrix(data)

# Create Groups
group = factor(c(
  rep("A_Wt", 10),
  rep("B_Tg", 13),
  rep("C_Proph_Ther_Rem", 3),
  rep("D_Ther_Rem", 10),
  rep("E_Ther_Hum", 10),
  rep("F_Ther_Enb", 10),
  rep("G_Ther_Cim", 10)
))


# Create empty dataframe
anova_table = data.frame()

# Recursive parse all genes
for( i in 1:length( matrixdata[ , 1 ] ) ) {
  
  # Create dataframe for each gene
  df = data.frame("gene_expression" = matrixdata[ i , ], 
                  "group" = group)
  
  # Apply ANOVA for gene i
  gene_aov = aov( gene_expression ~ group , data = df)
  
  # Apply tukey's post-hoc test on ANOVA results
  tukey = TukeyHSD( gene_aov , conf.level = 0.99)
  
  # Vector calling Tukey's values
  tukey_data = c(tukey$group["B_Tg-A_Wt", 1],
                 tukey$group["B_Tg-A_Wt", 4],
                 tukey$group["C_Proph_Ther_Rem-A_Wt", 1],
                 tukey$group["C_Proph_Ther_Rem-A_Wt", 4],
                 tukey$group["D_Ther_Rem-A_Wt", 1],
                 tukey$group["D_Ther_Rem-A_Wt", 4],
                 tukey$group["E_Ther_Hum-A_Wt", 1],
                 tukey$group["E_Ther_Hum-A_Wt", 4],
                 tukey$group["F_Ther_Enb-A_Wt", 1],
                 tukey$group["F_Ther_Enb-A_Wt", 4],
                 tukey$group["G_Ther_Cim-A_Wt", 1],
                 tukey$group["G_Ther_Cim-A_Wt", 4],
                 
                 tukey$group["C_Proph_Ther_Rem-B_Tg", 1],
                 tukey$group["C_Proph_Ther_Rem-B_Tg", 4],
                 tukey$group["D_Ther_Rem-B_Tg", 1],
                 tukey$group["D_Ther_Rem-B_Tg", 4],
                 tukey$group["E_Ther_Hum-B_Tg", 1],
                 tukey$group["E_Ther_Hum-B_Tg", 4],
                 tukey$group["F_Ther_Enb-B_Tg", 1],
                 tukey$group["F_Ther_Enb-B_Tg", 4],
                 tukey$group["G_Ther_Cim-B_Tg", 1],
                 tukey$group["G_Ther_Cim-B_Tg", 4])
  
  # Append Tukey's data to dataframe
  anova_table = rbind( anova_table , tukey_data)
}

colnames(anova_table) <- c("Wt_Tg_diff", "Wt_Tg_padj",
                           "Wt_Rem_P_diff", "Wt_Rem_P_padj", 
                           "Wt_Rem_diff", "Wt_Rem_padj", 
                           "Wt_Hum_diff", "Wt_Hum_padj", 
                           "Wt_Enb_diff", "Wt_Enb_padj", 
                           "Wt_Cim_diff", "Wt_Cim_padj",
                           
                           "Tg_Rem_P_diff", "Tg_Rem_P_padj", 
                           "Tg_Rem_diff", "Tg_Rem_padj", 
                           "Tg_Hum_diff", "Tg_Hum_padj", 
                           "Tg_Enb_diff", "Tg_Enb_padj", 
                           "Tg_Cim_diff", "Tg_Cim_padj")

# Add rownames with gene names
rownames(anova_table) = Gene

8 Volcano Plot

Code 
# -----------------Volcano plot preparation---------------------

# Create variables for upregulated/downregulated genes and genes with no observed change in expression levels
upWT = 0
downWT = 0
nochangeWT = 0

# Filter Differential Expressed Genes
upWT = which(anova_table[ , 1 ] < -1.0 & anova_table[ , 2 ] < 0.05)

downWT = which(anova_table[ , 1 ] > 1.0 & anova_table[ , 2 ] < 0.05)

nochangeWT = which(anova_table[ , 2 ] > 0.05 | 
                  (anova_table[ , 1 ] > -1.0 & anova_table[ , 1 ] < 1.0 ) )

# Create vector to store states for each gene
state <- vector(mode="character", length=length(anova_table[,1]))
state[upWT]   <- "up_WT"
state[downWT] <- "down_WT"
state[nochangeWT] <- "nochange_WT"

# Identify names of genes differentially expressed between wt and tg
genes_up_WT   <- c(rownames(anova_table)[upWT])
genes_down_WT <- c(rownames(anova_table)[downWT])

# Union of DEGs between wt and tg
deg_wt_tg <- c(genes_up_WT, genes_down_WT)

# Subset dataframe based on specific degs
deg_wt_tg_df <- subset(data , Gene %in% deg_wt_tg)

## Dataframe for volcano plot
volcano_data <- data.frame("padj" = anova_table[,2], 
                           "DisWt" = anova_table[,1], 
                           state=state)
Code 
ggplot(volcano_data , aes(x = DisWt , y = -log10(padj) , colour = state )) +
    
    geom_point() +
    
    labs(x = "mean(Difference)",
         y = "-log10(p-value)",
         title = "Volcano Plot",
         subtitle = "Differentially Expressed Genes (WT vs TG)") +
    
    # Insert line to show cutoff
    geom_vline(xintercept = c( -1 , 1 ),
               linetype = "dashed",
               color = "black") +
    
    # insert line to show cutoff
    geom_hline(yintercept = -log10(0.05),
               linetype = "dashed",
               color = "black")
Figure 8: Volcano plot of differentially expressed genes

8.1 Notes on Volcano Plot Interpretation

A volcano plot is a type of scatter plot commonly used in high-throughput data analysis, particularly in genomics, proteomics, and transcriptomics. It’s particularly useful for visualizing the results of differential expression analyses, where you’re comparing two conditions to identify genes, proteins, or other features that are significantly differentially expressed.

Structure

  • X-axis (mean(Difference)): This represents the magnitude of change between two conditions. A mean change of 0 indicates no change, positive values indicate upregulation (more abundant in the condition being tested), and negative values indicate downregulation (less abundant).

  • Y-axis (−log10(p-value)): This axis represents the statistical significance of the observed changes, often using a p-value from a statistical test. The higher the value on this axis, the more significant the change. Since it is a negative log scale, higher values represent smaller p-values.

Interpretation

  • Significance: Points located at the top of the plot have high statistical significance. The further a point is from the origin along the y-axis, the more significant the result.

  • Magnitude: Points farther to the right (positive mean change) indicate genes/proteins/features that are upregulated, while points farther to the left (negative mean change) indicate downregulation.

  • Significant Features: Typically, the most interesting features are those that are both statistically significant (high on the y-axis) and have a large magnitude of change (far left or far right on the x-axis). These appear as points in the upper left or upper right corners of the plot.

  • Thresholds: Horizontal and vertical lines are often added to represent thresholds for significance. Points outside these thresholds are often colored differently to highlight significant changes.

9 Uniform Manifold Approximation and Projection (UMAP) after identifying Differentially Expresssed Genes

Code 
# Subset dataframe based on specific degs
deg_wt_tg_df = subset(data , Gene %in% deg_wt_tg)

deg_wt_tg_df = deg_wt_tg_df[,1:23]

# After dataframe transposition columns must represent genes
deg_wt_tg_df = t(deg_wt_tg_df)
Code 
# UMAP dimension reduction for wt and tg samples
deg_wt_tg_df.umap = umap(deg_wt_tg_df, n_components=2, random_state=15)

# Keep the numeric dimensions
deg_wt_tg_df.umap = deg_wt_tg_df.umap[["layout"]]

# Create vector with groups
group = c(rep("A_Wt", 10), rep("B_Tg", 13) )

# Create final dataframe with dimensions and group for plotting
deg_wt_tg_df.umap = cbind(deg_wt_tg_df.umap, group)
deg_wt_tg_df.umap = data.frame(deg_wt_tg_df.umap)

# Plot UMAP results
ggplotly(
  ggplot(deg_wt_tg_df.umap, aes(x = V1, y = V2, color = group)) +
    geom_point() +
    labs(x = "UMAP1", y = "UMAP2", 
         title = "UMAP plot", 
         subtitle = "A UMAP Visualization of WT and TG samples (DEGs subset)") +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks = element_blank())
  )
Figure 9: UMAP plot of differentially expressed genes
Code 
# Group wt and tg as character and not factor
group = c(rep("A_Wt", 10), rep("B_Tg", 13) )

# Dimension reduction with PCA for wt and tg dataframe
deg_wt_tg_df.pca = prcomp(deg_wt_tg_df , scale. = FALSE)

deg_wt_tg_df.pca = data.frame("PC1" = deg_wt_tg_df.pca$x[,1] , 
                              "PC2" = deg_wt_tg_df.pca$x[,2] , 
                              "group" = group)

# Plot PCA results
ggplotly(
  ggplot(deg_wt_tg_df.pca , aes(x=PC1,y=PC2,color=group))+
    geom_point()+
    labs(x = "PC1", y = "PC2", 
         title = "PCA plot", 
         subtitle = "A PCA Visualization of WT and TG samples (DEGs subset)") +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks = element_blank())
)
Figure 10: PCA plot of differentially expresssed genes

9.0.1 Identify differentially expressed genes between transgenic animals and at least one therapy

Code 
# Volcano plot dataframe preparation for DEGs from TG vs therapies
upTHER = 0
downTHER = 0
nochangeTHER = 0

# Filter genes based on mean diff and p_value between TG and therapies
upTHER = which((anova_table[,13] < -1.0 & anova_table[,14] < 0.05) | 
               (anova_table[,15] < -1.0 & anova_table[,16] < 0.05) | 
               (anova_table[,17] < -1.0 & anova_table[,18] < 0.05) |
               (anova_table[,19] < -1.0 & anova_table[,20] < 0.05) |
               (anova_table[,21] < -1.0 & anova_table[,22] < 0.05) )

downTHER = which((anova_table[,13] > 1.0 & anova_table[,14] < 0.05) | 
                 (anova_table[,15] > 1.0 & anova_table[,16] < 0.05) | 
                 (anova_table[,17] > 1.0 & anova_table[,18] < 0.05) |
                 (anova_table[,19] > 1.0 & anova_table[,20] < 0.05) |
                 (anova_table[,21] > 1.0 & anova_table[,22] < 0.05) )

nochangeTHER = which( ( (anova_table[,13] > -1.0 & anova_table[,13] < 1.0) |   
                         anova_table[,14] > 0.05) |
                  
                      ( (anova_table[,15] > -1.0 & anova_table[,15] < 1.0) |
                         anova_table[,16] > 0.05) |
                        
                      ( (anova_table[,17] > -1.0 & anova_table[,17] < 1.0) |
                         anova_table[,18] > 0.05) |
                        
                      ( (anova_table[,19] > -1.0 & anova_table[,19] < 1.0) |
                         anova_table[,20] > 0.05) |
                        
                      ( (anova_table[,21] > -1.0 & anova_table[,21] < 1.0) |
                         anova_table[,22] > 0.05) )

# Create vector to store states for each gene
state = vector(mode = "character", length = length(anova_table[,1]))
state[upTHER] = "up_THER"
state[downTHER] = "down_THER"
state[nochangeTHER] = "nochange_THER"

# Identify names of genes differentially expressed between TG and therapies
genes_up_THER = c(rownames(anova_table)[upTHER])
genes_down_THER = c(rownames(anova_table)[downTHER])

deg_tg_ther = c(genes_up_THER, genes_down_THER)

# Combine DEGs from TG and THER
DEGs = c(deg_tg_ther, deg_wt_tg)

# Data frame with all DEGs for clustering
DEGsFrame = anova_table[rownames(anova_table) %in% DEGs, ]
DEGsFrame = as.matrix(DEGsFrame)

10 Hierarchical clustering

Code 
# k-means clustering
# ------------------
kmeans = kmeans(DEGsFrame[, c(1,3,5,7,9,11) ], centers = 6)
ggplotly(fviz_cluster(kmeans, data = (DEGsFrame[, c(1,3,5,7,9,11)]), geom = "point", show.clust.cent = TRUE))
Figure 11: Clusters of differentially expressed genes
Code 
# Extract genes from clusters
clusters = data.frame(kmeans$cluster)
colnames(clusters) = ("ClusterNo")

# Output data as table (group by cluster number)
kable(clusters, col.names = "Cluster Number", caption = "Clusters and associated genes") |>
        kable_styling(font_size = 14) |>
        scroll_box(height = "400px")
Clusters and associated genes
Cluster Number
Abca5 1
Abca6 1
Abca8a 1
Abca8b 1
Abcb4 4
Abhd14b 5
Ablim1 1
Abra 2
Acaca 1
Acacb 4
Acan 1
Ache 2
Acp2 3
Acp5 6
Acpp 6
Acsl6 4
Acss3 1
Acta1 2
Actc1 1
Actg2 5
Actn2 2
Actn3 2
Acvr1c 1
Acvr2a 5
Acyp2 2
Adam8 6
Adamts13 3
Adamts15 3
Adamts4 6
Adamts5 1
Adamts7 3
Adcy5 5
Adig 1
Adipoq 1
Adk 1
Adprhl1 2
Adrbk2 6
Aff2 4
Agbl1 2
Agpat9 5
Aif1 3
Aim1 6
Akna 3
Akt2 1
Alas2 4
Alcam 6
Aldh1a1 1
Aldh1a3 1
Aldh1a7 1
Aldh3a1 5
Aldh6a1 1
Aldob 3
Aldoc 3
Alg3 3
Alox12 1
Alpk3 2
Amacr 3
Amot 2
Ampd1 2
Amy1 1
Angptl4 3
Angptl7 1
Ank2 5
Ank3 2
Ankrd1 2
Ankrd2 4
Ankrd23 2
Ankrd55 3
Ano5 4
Anpep 3
Anxa8 5
Aoah 6
Aox1 1
Aox3 5
Ap4s1 1
Apbb1ip 3
Apobec1 6
Apobec2 4
Apobec3 3
Aprt 3
Aqp7 1
Ar 2
Arg1 6
Arg2 3
Arhgap19 1
Arhgap20 4
Arhgap28 5
Arhgap30 6
Arhgap4 3
Arhgef10 5
Arhgef9 1
Arid3c 3
Arpc1b 3
Arpc2 3
Arpc3 3
Arrdc4 3
Art3 4
Asb11 2
Asb12 4
Asb14 4
Asb15 2
Asb2 2
Asb5 2
Asf1b 3
Aspa 4
Aspm 3
Asrgl1 5
Atf3 3
Atn1 3
Atp1a2 4
Atp1a3 3
Atp1b1 2
Atp1b2 2
Atp1b4 2
Atp2a1 1
Atp6v0a4 1
Atp6v0d2 6
Atp6v1a 3
Atp6v1b1 3
Atp8b4 6
Atp9a 1
Atpbd4 3
Aurkb 3
B3galt2 4
B3gnt5 3
B4galnt1 3
B4galt5 3
B4galt6 3
Baiap2 6
Bak1 3
Bank1 3
Basp1 3
Bax 6
Bbox1 4
Bccip 5
Bche 4
Bcl10 3
Bcl11b 3
Bcl2a1a 3
Bcl2a1b 3
Bcl2a1d 3
Bcl2l1 5
Bcl3 6
Bdh2 1
Bdkrb1 6
Bend4 3
Best3 2
Bid 3
Bin1 2
Birc3 6
Blk 3
Blnk 3
Blvrb 3
Bmper 5
Bod1 1
Brms1l 5
Bst1 6
Btg1 5
Btg2 3
Btk 6
Btla 3
Bub1 3
Bub1b 6
Bves 4
C1qb 3
C1qtnf3 6
C1rl 3
C3 6
C4a 5
C4b 1
C7 4
C77370 1
C87977 5
Cacna1i 3
Cacna1s 2
Cacna2d1 2
Cacnb1 2
Cacng1 2
Cadm3 2
Calcr 3
Camk2a 2
Camk2b 4
Camkk1 4
Cand2 2
Cap2 4
Capn3 4
Capn6 3
Car12 6
Car3 1
Car5b 5
Car8 1
Casc5 3
Casp1 6
Casp4 6
Casp7 3
Casp8 3
Casq1 2
Casq2 1
Casr 1
Cav3 4
Cbr2 4
Ccbp2 1
Ccdc23 3
Ccdc3 4
Ccdc62 5
Ccdc67 3
Ccl11 2
Ccl17 6
Ccl19 6
Ccl2 6
Ccl20 3
Ccl21a 6
Ccl22 6
Ccl24 4
Ccl28 3
Ccl3 3
Ccl5 6
Ccl7 6
Ccl8 3
Ccl9 6
Ccna2 3
Ccnb1 6
Ccnb2 3
Ccnd1 3
Ccne2 3
Ccr2 6
Ccr3 3
Ccr6 3
Ccr7 3
Ccr8 3
Ccrl2 3
Cd101 3
Cd14 6
Cd163 1
Cd177 5
Cd19 3
Cd22 3
Cd248 5
Cd27 3
Cd274 3
Cd300lb 3
Cd33 3
Cd37 6
Cd38 3
Cd3d 3
Cd3e 3
Cd3g 6
Cd4 6
Cd44 6
Cd46 5
Cd48 3
Cd5 3
Cd52 6
Cd53 6
Cd55 1
Cd59a 2
Cd63 3
Cd68 6
Cd69 3
Cd72 6
Cd74 3
Cd79a 6
Cd79b 3
Cd80 6
Cd84 6
Cdadc1 5
Cdc25c 3
Cdh19 4
Cdh22 3
Cdh5 3
Cdhr1 3
Cdk18 3
Cdk2ap2 3
Cdk6 6
Cdkl1 3
Cdkn1c 5
Cdo1 5
Cdon 1
Ceacam1 6
Ceacam16 3
Ceacam19 6
Ceacam2 3
Cenpe 6
Cenpf 3
Cenpm 3
Cenpn 3
Cep55 3
Cetn4 5
Cfb 6
Cfd 1
Cfp 3
Ch25h 5
Chad 4
Chd2 5
Chi3l1 6
Chl1 6
Chodl 4
Chpt1 4
Chrna1 2
Chrnb1 2
Chst2 1
Cidec 4
Ciita 6
Cilp 5
Cilp2 2
Cirbp 1
Ckap2 3
Ckap2l 3
Ckb 3
Ckmt2 4
Clcn1 2
Cldn20 5
Clec12a 6
Clec2d 3
Clec3a 4
Clec3b 2
Clec4a1 6
Clec4a2 3
Clec4a3 6
Clec4d 6
Clec4e 6
Clec4n 6
Clec5a 6
Clec7a 6
Clic5 4
Cln5 3
Cln6 3
Clnk 3
Clu 5
Cmbl 4
Cmya5 2
Cnksr1 2
Cnn2 3
Cnn3 3
Cnr2 3
Cobll1 5
Coch 4
Col10a1 4
Col12a1 6
Col22a1 5
Col28a1 5
Comp 1
Coro1a 3
Coro2a 3
Coro6 2
Cotl1 3
Cox6a2 2
Cox7a1 2
Cox8a 3
Cox8b 2
Cpa3 2
Cpa6 3
Cpe 5
Cpeb3 4
Cpm 4
Cpne8 5
Cpxm1 6
Cr2 3
Crb3 3
Creb5 2
Crebl2 5
Creg1 3
Crispld1 1
Cry2 1
Cryaa 3
Cryab 4
Cryba2 3
Csf1r 3
Csf2ra 6
Csf2rb 6
Csf2rb2 6
Csf3r 3
Csgalnact2 5
Csmd1 6
Csrp3 2
Cst7 3
Cstb 3
Cth 1
Cthrc1 6
Ctla2b 3
Ctnnal1 1
Cts3 5
Ctsc 3
Ctsk 6
Ctss 6
Ctsz 6
Cttnbp2 5
Ctxn1 3
Ctxn3 4
Cubn 3
Cx3cl1 6
Cx3cr1 5
Cxcl1 6
Cxcl10 3
Cxcl13 6
Cxcl16 6
Cxcl2 6
Cxcl3 6
Cxcl5 6
Cxcl9 3
Cxcr2 3
Cxcr4 6
Cyba 3
Cybb 6
Cyp2e1 1
Cyp2f2 1
Cyp2s1 6
Cyp3a13 1
Cyp4f18 3
Cyp4v3 6
Cyp7b1 6
Cypt4 3
Cyth4 6
Cytip 6
Cytl1 4
Daam2 1
Dad1 3
Dapp1 3
Dcbld2 5
Ddit4l 4
Ddrgk1 5
Ddx60 5
Des 2
Dgat2 4
Dgkb 1
Dhrs7c 1
Diap3 6
Dixdc1 1
Dkk2 4
Dmd 2
Dmp1 6
Dmpk 4
Dmxl2 6
Dna2 3
Dnaja4 4
Dnajc21 5
Dnajc22 3
Dock10 6
Dock2 6
Dock5 6
Dok3 6
Dok6 3
Dpep2 6
Dpf3 2
Dpm1 5
Dram1 6
Drp2 4
Dsg1b 3
Dsg2 3
Dstyk 5
Dtl 3
Dtna 2
Dtx3l 3
Dusp10 2
Dusp13 2
Dusp14 1
Dusp27 2
Dusp5 3
Dync2h1 5
E2f1 3
E2f4 3
E2f8 3
Ear2 3
Ebf2 1
Ebi3 3
Ebna1bp2 5
Ecscr 3
Ect2 6
Eda2r 4
Edil3 6
Ednrb 2
Eef1a2 4
Efcab2 4
Efha2 1
Efhd1 1
Efhd2 6
Egf 5
Egln3 2
Egr2 3
Ehbp1 1
Ehd4 5
Eif3k 3
Eln 5
Eme1 3
Emilin1 6
Emr1 6
Emr4 6
Emx2 3
Eng 3
Eno3 2
Enpep 1
Entpd7 3
Epb4.1l4a 2
Epb4.1l4b 1
Epha3 3
Epm2a 2
Eps8l2 5
Epsti1 3
Espl1 3
Etv6 3
Evc 5
Evi2a 3
Evi5 5
Exo1 3
Exosc1 3
Ext1 6
Eya4 2
F10 6
F13a1 1
F2r 3
F2rl1 3
F2rl2 3
F7 6
Fabp4 1
Fabp7 6
Fabp9 5
Fads3 3
Faim3 3
Fam105a 6
Fam108c 3
Fam114a2 5
Fam124b 3
Fam126a 5
Fam126b 1
Fam129a 1
Fam13a 4
Fam160a1 4
Fam35a 5
Fam64a 3
Fam83f 3
Fancc 5
Fas 6
Fasn 1
Fbl 3
Fbn2 5
Fbp2 3
Fbxo30 5
Fbxo40 4
Fcer1g 6
Fcer2a 3
Fcf1 5
Fcgr2b 6
Fcgr3 6
Fcgr4 3
Fcrl1 3
Fcrla 3
Fcrlb 3
Fermt2 1
Fermt3 6
Fez1 1
Fgf13 4
Fgf2 1
Fgf23 3
Fgf9 1
Fgfr3 5
Fgl2 1
Fgr 6
Fhl1 1
Fhl3 2
Fhod3 2
Figf 5
Filip1 4
Fitm1 4
Fkbp2 5
Flnc 2
Flt3 6
Fmnl1 3
Fmo2 1
Fndc5 2
Folr2 2
Fos 3
Fosl2 3
Foxp3 3
Fpr1 3
Fras1 2
Frmd3 1
Frmd7 1
Frrs1 6
Frzb 1
Fscn1 6
Fsd1l 4
Fsd2 4
Fv1 5
Fxyd1 2
Fxyd2 3
Fyb 6
Fyco1 2
G0s2 1
Gabra4 1
Gabrb3 3
Galnt6 6
Galntl2 1
Gbe1 4
Gbgt1 6
Gbp2 6
Gbp5 3
Gclm 3
Gdap1 2
Gdf10 1
Gdpd1 3
Gfpt2 1
Gfra2 3
Gimap4 5
Gimap5 3
Gjb3 6
Gjc3 2
Gk5 5
Gla 6
Glb1l2 1
Gldn 4
Gli1 1
Glrx 3
Gls 5
Glt25d2 4
Glt28d2 4
Glycam1 3
Gm10228 6
Gm10229 6
Gm11559 3
Gm14492 2
Gm1679 3
Gm2a 3
Gm3893 5
Gm4792 3
Gm5150 6
Gm527 5
Gm5465 3
Gm5563 3
Gm6026 3
Gm614 6
Gm7609 6
Gm885 6
Gm889 2
Gm9733 3
Gm9766 1
Gmpr 2
Gna15 3
Gnai1 1
Gngt2 3
Gnpnat1 1
Gnptab 6
Gp49a 6
Gpam 2
Gpc3 4
Gpc4 1
Gpd1 1
Gphb5 3
Gpr1 4
Gpr114 3
Gpr132 6
Gpr137b 6
Gpr141 3
Gpr161 5
Gpr162 3
Gpr165 4
Gpr176 6
Gpr18 3
Gpr35 3
Gpr64 1
Gpr65 6
Gpr68 3
Gpr84 6
Gpr97 3
Gprc5a 2
Gpt2 2
Grap2 3
Grb14 1
Gria2 5
Grid2 1
Grina 3
Gsn 5
Gsta3 5
Gsta4 1
Gstk1 2
Gusb 6
Gxylt2 6
Gyg 2
Gypa 3
Gys1 2
Gzmc 3
H19 1
Hacl1 5
Hal 3
Hapln1 5
Has1 1
Has2 1
Havcr2 6
Hbegf 1
Hck 6
Hcls1 6
Hcn1 1
Hectd2 1
Hemgn 3
Hexa 3
Hfe2 2
Hhatl 2
Hhip 4
Hibch 5
Higd1a 5
Hist1h2ab 6
Hist1h2bb 3
Hist2h2bb 3
Hmcn1 6
Hmgb2 5
Hmha1 3
Hmox1 6
Hn1 3
Hn1l 3
Hnrnpm 5
Homer2 2
Hoxa11 5
Hoxa13 2
Hoxc10 1
Hoxd13 2
Hp 6
Hpse2 2
Hrc 2
Hspa12a 1
Hspb1 4
Hspb3 4
Hspb6 4
Hspb7 2
Hspb8 4
Htatsf1 5
Htra2 3
Htra4 1
Hvcn1 6
Icam1 6
Id4 1
Idua 3
Ifi205 3
Ifi27l2a 1
Ifi30 6
Ifitm3 3
Ifitm6 1
Ifnar1 6
Ifnar2 3
Ifngr2 3
Igf2 1
Igfbp3 3
Igfbp6 1
Igfbp7 3
Igj 3
Igsf10 6
Igsf6 6
Ikbke 6
Ikzf3 3
Ikzf5 5
Il10ra 6
Il12rb2 3
Il13ra1 6
Il16 5
Il1a 6
Il1b 6
Il1f9 3
Il1rn 6
Il20ra 3
Il23r 3
Il28ra 3
Il2ra 3
Il2rg 6
Il6 6
Il7r 6
Inmt 4
Insig2 5
Insl6 3
Ipp 5
Iqgap3 3
Irak2 3
Irg1 6
Itga2 3
Itga4 6
Itga5 6
Itga6 1
Itga7 4
Itga8 1
Itgae 3
Itgal 6
Itgam 6
Itgax 6
Itgb1bp2 2
Itgb2 6
Itgb3 6
Itgb6 4
Itgb8 1
Itih2 1
Itih5 5
Jdp2 6
Jph1 2
Jph2 2
Jup 2
Kank1 5
Kbtbd10 4
Kbtbd5 2
Kcna1 2
Kcna3 6
Kcna6 5
Kcnc1 4
Kcnc2 1
Kcne1 3
Kcng3 3
Kcnj2 4
Kcnn4 6
Kcnq5 1
Kcnt2 1
Kctd2 5
Kdm5a 5
Kif11 3
Kif13b 5
Kif1b 5
Kif20a 6
Kif21a 1
Kif26b 3
Klf12 1
Klf4 2
Klhl13 1
Klhl31 2
Klra17 3
Klra2 3
Klrb1f 3
Klrc1 3
Klrd1 3
Kmo 3
Kng2 6
Krt31 3
Krt6a 3
Krtap1-5 6
Krtap3-2 3
Krtap4-1 3
Krtdap 4
Ky 2
Kynu 3
Lair1 6
Lama2 4
Lamb2 1
Lass6 6
Lat2 6
Layn 3
Lcp1 6
Lcp2 6
Ldb3 2
Lef1 3
Lep 4
Lgals12 4
Lgi1 4
Lgi4 5
Lif 6
Lilra6 6
Lilrb3 3
Lilrb4 6
Limch1 2
Lims2 4
Lipa 6
Lipn 6
Lmcd1 2
Lmo2 5
Lmod2 2
Lmod3 4
Lpar1 5
Lpcat2 6
Lpin1 2
Lpl 1
Lrdd 3
Lrrc15 6
Lrrc2 4
Lrrc25 3
Lrrc39 4
Lrrc55 3
Lrrc61 3
Lrrn1 4
Lrrn4cl 2
Lrtm1 4
Ltb 6
Ltbp2 6
Luc7l3 5
Lxn 6
Ly6g5c 3
Ly6i 5
Ly86 3
Ly9 6
Ly96 3
Lyn 3
Lynx1 2
Lyve1 1
Lyz1 6
Lyz2 6
Maff 3
Mageh1 5
Malt1 6
Mamdc2 5
Manea 3
Map1lc3a 4
Map4k1 3
Mapk11 3
Mapk12 4
Mapk3 3
Mapkapk3 3
Mapt 2
Marcksl1 5
Marco 6
Matn2 5
Mb 2
Mccc1 5
Mcm5 3
Mcoln2 6
Mcpt4 4
Mdga2 1
Med29 3
Mef2c 1
Mefv 6
Meg3 1
Megf10 1
Meox1 3
Mfap4 3
Mfi2 1
Mfn2 1
Mfsd4 3
Mgll 1
Mgst3 2
Micall2 3
Mif 3
Mir100 1
Mir125b-1 1
Mir133a-1 2
Mir133a-2 2
Mir133b 4
Mir142 6
Mir194-2 5
Mir23b 5
Mir29b-2 1
Mir300 5
Mir376b 5
Mir380 1
Mir382 5
Mir487b 5
Mir543 5
Mirlet7a-2 5
Mirlet7c-1 1
Mirlet7c-2 5
Mki67 6
Mlf1 2
Mlxip 3
Mlxipl 2
Mme 3
Mmp12 6
Mmp13 6
Mmp14 6
Mmp19 6
Mmp25 3
Mmp3 6
Mmp9 6
Mn1 4
Mobp 3
Mocs1 3
Mpdz 5
Mpeg1 6
Mpp3 4
Mpp7 5
Mpz 4
Mpzl3 3
Mreg 2
Mrgprb1 4
Mrgprb2 4
Mrpl22 3
Mrpl34 3
Mrpl54 3
Mrps24 3
Ms4a1 3
Ms4a14 6
Ms4a6d 6
Ms4a7 6
Msc 3
Msr1 6
Mstn 2
Mt2 3
Mtap1b 2
Mtmr11 5
Murc 2
Musk 4
Mustn1 2
Myadml2 4
Mybpc1 2
Mybpc2 2
Mybph 2
Myc 3
Myf6 4
Myh1 1
Myh10 1
Myh11 5
Myh2 4
Myh3 1
Myh4 5
Myh7 2
Myh8 1
Myl1 2
Myl2 4
Myl3 2
Myl6b 2
Mylk2 2
Mylk4 4
Mylpf 2
Myo18b 4
Myo1f 6
Myom1 1
Myom2 1
Myom3 4
Myot 2
Myoz1 2
Myoz2 2
Myoz3 4
Mypn 2
N4bp1 3
Nagpa 3
Naip2 6
Naip5 6
Naip6 6
Napsa 6
Nat14 5
Ncf1 6
Ncf2 6
Ncf4 6
Nckap1l 6
Ncr1 3
Ndc80 6
Ndrg2 4
Ndufb9 3
Ndufs8 2
Neb 2
Necab1 4
Neil3 3
Nes 4
Neurl3 6
Nexn 4
Nfam1 6
Nfatc1 3
Nfe2 3
Nfe2l2 3
Nfkb2 6
Nfkbia 6
Nfkbid 6
Nfkbie 6
Ngf 3
Nhedc2 6
Ninl 3
Nkg7 3
Nlrc4 3
Nlrp3 6
Nmral1 3
Nnat 4
Nop10 3
Nos2 6
Nov 1
Nova1 4
Npr1 1
Npr2 1
Npr3 4
Npy 4
Nr1d1 4
Nr2f2 5
Nr3c2 1
Nr4a3 3
Nr5a2 3
Nrap 2
Nrbf2 3
Nrg4 1
Nrn1 5
Nrn1l 4
Nrp2 6
Nt5dc2 6
Ntn1 5
Ntn4 1
Ntng1 1
Ntrk2 4
Ntrk3 4
Nudt10 5
Nuf2 6
Nupl1 3
Nxph1 3
Obp1a 5
Obscn 1
Ociad2 2
Olfr1020 5
Olfr1130 5
Olfr1249 5
Olfr1252 5
Olfr1307 5
Olfr133 5
Olfr1413 5
Olfr1474 5
Olfr172 5
Olfr179 3
Olfr444 5
Olfr605 3
Olfr684 5
Olfr732 5
Olfr781 5
Olfr866 5
Olfr888 5
Olfr889 5
Olfr921 5
Olfr963 2
Olr1 6
Omd 2
Osbpl1a 5
Osbpl6 2
Oscar 6
Osm 3
Ostn 4
Otud1 2
Ovgp1 5
Oxsr1 5
P2rx1 5
P2rx4 3
P2ry10 3
P2ry6 3
P4ha3 6
Pacsin3 1
Palmd 2
Panx1 3
Panx3 6
Pcbp4 2
Pcdh11x 4
Pcdh9 1
Pck1 1
Pcolce2 4
Pcp4l1 2
Pcsk6 1
Pcx 1
Pdcd1lg2 3
Pde1a 5
Pde4dip 2
Pdha1 5
Pdk2 2
Pdk4 3
Pdlim3 2
Pdlim4 3
Pdpn 6
Peg3 1
Penk 4
Pfkfb1 4
Pfkfb3 5
Pfkm 4
Pfn2 4
Pgam2 4
Pgm2 1
Pgm3 5
Phka1 2
Phkg1 2
Phtf2 1
Pi15 1
Pi16 3
Pid1 3
Pik3ap1 6
Pik3cg 6
Pik3r5 6
Pilra 3
Pilrb1 3
Pion 6
Pkia 4
Pla1a 1
Pla2g16 1
Pla2g2a 4
Pla2g2d 6
Pla2g4e 2
Pla2g7 6
Plac8 3
Plaur 6
Plbd1 6
Plcb2 6
Plcg2 3
Pld3 6
Pld4 6
Plek 6
Plek2 6
Plin1 1
Plin4 2
Plk1 3
Plk3 3
Plp1 4
Pls3 5
Plscr1 5
Plxnb2 6
Plxnc1 3
Pmp2 4
Pmp22 5
Pnn 5
Pnpla3 4
Podn 2
Podnl1 6
Pold4 3
Popdc2 4
Popdc3 2
Postn 6
Pot1b 3
Pou2af1 3
Pou2f2 3
Pparg 1
Ppargc1a 2
Ppfibp2 5
Ppic 6
Ppip5k1 4
Ppl 1
Ppp1r14c 1
Ppp1r3a 2
Ppp1r3c 4
Ppp2r3a 2
Ppp4r4 5
Prc1 6
Prcp 6
Prdm1 6
Prdm13 3
Prelid1 3
Prkaa2 4
Prkar2a 2
Prkcd 6
Prkcq 2
Prkd1 1
Prkg1 1
Prl2c5 3
Prnp 5
Prokr1 3
Prps1 5
Prr11 3
Prr23a 3
Prss46 3
Prss50 3
Prune2 2
Psd4 3
Psmb10 3
Psmb8 3
Psmd10 5
Pstpip1 6
Ptafr 6
Ptbp1 3
Ptch1 1
Ptger2 3
Ptger3 5
Ptger4 6
Ptgs1 1
Ptgs2 3
Ptk2b 6
Ptpla 2
Ptplad2 3
Ptpn22 6
Ptpn3 4
Ptpn6 6
Ptpn7 3
Ptprc 6
Ptpre 6
Ptprv 6
Ptrf 1
Ptx4 2
Pvalb 4
Pvrl1 6
Pycard 3
Pygm 2
Qpct 2
Rab11fip1 3
Rab32 6
Rab37 1
Rab38 6
Rabgap1l 5
Rac2 6
Racgap1 6
Rai14 6
Rasd1 1
Rasgrp1 6
Rassf2 3
Rassf4 6
Rassf6 5
Raver2 1
Rbm24 2
Rbm7 3
Rbp1 6
Rcan2 4
Rcn2 5
Reck 1
Rel 3
Relb 6
Rep15 5
Retn 2
Retnla 4
Rgs1 6
Rgs19 6
Rgs5 1
Rhag 3
Rhbdf2 3
Rhog 3
Rilpl2 3
Rin3 3
Rinl 3
Rnaset2a 3
Rnf122 3
Rnf125 3
Rnf128 6
Rnf146 5
Rnf149 6
Rnf157 3
Rnf19b 6
Rnu3a 3
Rnu73b 3
Rny1 6
Rora 5
Rorc 4
Rp2h 3
Rpa3 3
Rpl28 3
Rpl34 3
Rpl3l 2
Rps19 3
Rragb 1
Rragd 2
Rrs1 5
Rsu1 3
Rtn2 4
Rtn4rl1 1
Ruvbl1 5
Rxrg 4
Ryr1 4
S100a8 3
S100b 1
Saa3 6
Samsn1 3
Sash3 3
Satb1 1
Sbk2 4
Sbsn 1
Scd1 1
Scg3 4
Scg5 5
Scn4a 2
Scn4b 4
Scrg1 1
Sdf2l1 3
Sectm1b 5
Sel1l3 2
Sele 6
Sell 3
Selp 6
Selplg 6
Sema3a 1
Sema3c 2
Sema3d 1
Sema3e 1
Sema4a 6
Sema4d 3
Sema6c 4
Serpina3c 4
Serpina3f 6
Serpinb9 3
Serpine1 6
Serpine2 5
Sestd1 5
Sf3b5 3
Sfpi1 3
Sfrp4 5
Sfrp5 4
Sgca 2
Sgcd 1
Sgcg 2
Sgsh 3
Sgsm2 5
Sh2d1b1 3
Sh3bgrl3 3
Sh3bp1 3
Sh3bp2 3
Sh3kbp1 3
Sh3rf2 2
Siglece 3
Siglecg 3
Sipa1 3
Sirpb1a 6
Sirpb1b 6
Sirt6 3
Sla 3
Slain1 3
Slamf1 3
Slamf7 3
Slamf8 6
Slc10a6 6
Slc11a1 6
Slc13a3 3
Slc14a1 3
Slc15a3 6
Slc16a10 3
Slc16a6 3
Slc18a1 3
Slc19a3 3
Slc20a1 5
Slc22a3 1
Slc25a37 3
Slc25a4 1
Slc25a5 3
Slc2a4 2
Slc2a6 6
Slc31a2 3
Slc33a1 5
Slc36a1 3
Slc37a2 6
Slc38a3 2
Slc38a4 4
Slc39a4 3
Slc47a1 4
Slc48a1 3
Slc7a2 6
Slc7a9 5
Slc8a1 3
Slco2b1 5
Slfn2 6
Slfn9 3
Slpi 6
Slurp1 2
Smarca1 5
Smpdl3b 6
Smpx 4
Smtn 5
Smtnl1 2
Smtnl2 4
Smyd1 4
Snai2 3
Snap25 2
Sncg 1
Sned1 5
Snora23 3
Snora44 3
Snora74a 3
Snord104 3
Snord118 3
Snord57 3
Snord82 3
Snrpn 2
Snx10 6
Snx13 5
Snx18 3
Snx20 3
Soat1 3
Socs5 5
Sod3 3
Sorbs1 2
Sox5 5
Sox6 2
Sox9 2
Spag5 3
Sparcl1 1
Speer6-ps1 5
Spg20 5
Spib 3
Spn 6
Spna1 3
Spnb1 2
Spock2 4
Spp1 6
Sprr1b 3
Sprr2e 3
Src 3
Srfbp1 5
Srl 4
Srsy 5
Ssty1 5
Ssty2 5
St18 3
St3gal5 1
St3gal6 1
St6galnac4 3
St8sia4 3
St8sia6 6
Stac3 2
Stap1 3
Stc1 3
Stc2 5
Stk32b 1
Stk38l 5
Stox1 3
Stra6 3
Stx4a 5
Stxbp2 3
Sucla2 1
Sult1e1 1
Susd2 6
Susd3 3
Sv2b 2
Sv2c 1
Syde2 5
Syn1 3
Syne2 5
Synpo2l 2
Sypl2 4
Tacc2 2
Taf10 3
Tagap 6
Tagln 5
Tarm1 3
Tarsl2 4
Tbc1d20 3
Tbc1d2b 3
Tbc1d8b 5
Tc2n 6
Tcap 2
Tcea3 2
Tchh 5
Tcirg1 3
Tctex1d2 3
Tenc1 2
Tff1 3
Tgfbr3 1
Thbd 1
Thrb 1
Thrsp 1
Tiaf2 5
Tifab 3
Tigd4 2
Timp1 6
Timp4 4
Tk1 3
Tlcd1 2
Tln2 5
Tlr1 3
Tlr13 6
Tlr2 6
Tlr8 6
Tlr9 3
Tm4sf19 6
Tm6sf1 3
Tmeff2 2
Tmem100 1
Tmem117 4
Tmem121 3
Tmem123 3
Tmem134 3
Tmem141 3
Tmem173 6
Tmem176a 6
Tmem176b 3
Tmem179b 3
Tmem182 2
Tmem196 1
Tmem45b 4
Tmem56 2
Tmem97 3
Tmod4 2
Tnf 6
Tnfaip2 6
Tnfaip3 6
Tnfaip6 6
Tnfrsf11a 6
Tnfrsf11b 1
Tnfrsf13c 3
Tnfrsf14 6
Tnfrsf1b 6
Tnfrsf26 3
Tnfrsf9 6
Tnfsf10 3
Tnfsf15 3
Tnfsf9 3
Tnik 4
Tnip1 6
Tnip3 6
Tnn 6
Tnnc1 4
Tnnc2 2
Tnni1 2
Tnni2 4
Tnnt1 4
Tnnt2 2
Tnnt3 2
Tns3 3
Tnxb 1
Tomm6 3
Tor2a 3
Tpm1 5
Tpm2 2
Tpsab1 4
Tpsb2 4
Tpx2 6
Traf1 6
Traf3ip3 3
Tram2 3
Trappc2l 3
Trdn 2
Trem1 6
Trem2 6
Trem3 3
Treml4 6
Trf 1
Trhde 1
Trim54 1
Trim55 2
Trim63 2
Trim65 3
Trim72 4
Trpc3 5
Trpc4 3
Tshr 4
Tshz2 1
Tspan15 1
Tspan17 3
Tspan2 5
Tspan3 5
Tspan6 3
Tspan8 2
Ttc9 4
Ttll7 2
Ttyh3 3
Tuft1 5
Txlnb 2
Txnrd3 1
Tyms 3
Tyrobp 3
Uaca 1
Uap1 5
Ube2c 3
Ube2q1 3
Ube2ql1 3
Ucma 4
Ucp2 6
Ucp3 3
Ufsp1 2
Ugcg 3
Ugdh 5
Ugp2 2
Ugt8a 1
Unc13d 3
Unc93b1 6
Upp1 3
Usmg5 2
Usp13 4
Usp2 2
Usp53 1
Utrn 5
Uxs1 5
Vash2 3
Vasp 6
Vat1l 1
Vav1 6
Vcam1 6
Vcan 5
Vgll2 2
Vgll3 3
Vit 1
Vkorc1 3
Vldlr 1
Vmn1r219 5
Vmn1r220 5
Vmn2r15 5
Vmn2r75 5
Vps37a 2
Vsig4 4
Vwa1 5
Was 3
Wbscr25 3
Wdfy4 6
Wfdc1 2
Wif1 5
Wipf2 3
Wisp1 6
Xcl1 3
Xirp1 2
Xirp2 1
Xkr7 3
Xpr1 3
Yipf7 2
Zadh2 5
Zbp1 6
Zbtb10 5
Zbtb16 2
Zbtb41 5
Zc3h12a 3
Zfp106 1
Zfp212 5
Zfp330 5
Zfp36 3
Zfp365 5
Zfp36l2 3
Zfp385b 1
Zfp511 3
Zfp628 3
Zfp800 3
Zfp870 5
Zfpm2 3
Zfyve9 5
Zmiz2 3
Zmynd15 6
Zmynd17 4
Znhit6 5
Code 
# Extract clusters
cluster1 = rownames(subset(clusters, ClusterNo == 1))
cluster2 = rownames(subset(clusters, ClusterNo == 2))
cluster3 = rownames(subset(clusters, ClusterNo == 3))
cluster4 = rownames(subset(clusters, ClusterNo == 4))
cluster5 = rownames(subset(clusters, ClusterNo == 5))
cluster6 = rownames(subset(clusters, ClusterNo == 6))
Code 
# ---------------------Prepare data for heatmap-------------------------

# Define custom function to perform hierarchical clustering with the Ward.D2 linkage method
hclustfunc = function(x)
  hclust(x, method = "ward.D2")

# Define custom function to calculate pairwise Euclidean distances between data points
distfunc = function(x)
  dist(x, method = "euclidean")

# Perform clustering on rows and columns 
cl.row = hclustfunc(distfunc(DEGsFrame[, c(1,3,5,7,9,11)]))

# Extract cluster assignments of rows
gr.row = cutree(cl.row, k=6)

# Apply a set of color palette
colors = brewer.pal(5, "Set3")

heatmap <- heatmap.2(
          DEGsFrame[, c(1,3,5,7,9,11)],
          col = bluered(100), # blue-red color palette
          tracecol="black",
          density.info = "none",
          labCol = c("TG", "REM_P", "REM", "HUM", "ENB","CIM"),
          scale="none", 
          labRow="", 
          vline = 0,
          mar=c(6,2),
          RowSideColors = colors[gr.row],
          hclustfun = function(x) hclust(x, method = 'ward.D2')
)
Figure 12: Heatmap of differentially expressed genes after hierarchical clustering

10.1 Notes on Heatmap Interpretation

Rows: Each row represents a gene and the genes are clustered hierarchically (dendrogram on the left). Genes with similar expression patterns are grouped together.

Columns: They represent experimental conditions/ sample types.

Color Gradient (Color Key):

  • Red represents upregulated genes (positive values).

  • Blue represents downregulated genes (negative values).

  • The scale ranges from -4 to +4, indicating the intensity of differential expression (e.g. mean change).

Clustering:

  • The dendrogram on the top shows how the conditions are related based on the gene expression patterns. Conditions that cluster closer together (e.g., REM and TG) are likely to have similar expression profiles.

  • The dendrogram on the left shows how genes are grouped based on similarity in expression across the conditions.

Key Differences of Hierarchical vs K-means Clustering
Aspect Hierarchical Clustering K-means Clustering
Approach Builds hierarchy (Dendrogram) Partitions into flat clusters
Number of Clusters Determined post-analysis Must be predefined
Cluster Shape Can handle various shapes Assumes spherical clusters
Distance Metric

Multiple metrics

(e.g. Euclidean, Manhattan)

 Only Euclidean distance
Scalability Not scalable for large datasets Efficient and scalable
Visual Output Dendrogram Cluster assignment
Handling Outliers Sensitive to outliers Sensitive, but manageable

Performing k-means clustering followed by creating a heatmap with hierarchical clustering combines the strengths of both clustering methods to get a more insightful and comprehensive visualization of the data. This hybrid approach can be particularly useful when working with large and complex datasets, such as gene expression data. This combination leverages the strengths of both methods: k-means for global structure and hierarchical clustering for local, detailed exploration.

11 Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)

Once DEGs are identified, they need to be annotated to determine their biological functions, cellular localization, molecular interactions, and involvement in various biological pathways. This is often achieved by comparing DEGs to databases of known gene annotations, such as Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG).

Pathway analysis focuses on identifying interconnected networks of genes that collaborate to carry out specific biological functions or participate in common signaling pathways. This involves mapping DEGs onto existing biological pathways and identifying key regulatory nodes or hub genes within these pathways. Pathway analysis provides insights into the underlying molecular mechanisms driving the observed gene expression changes.

The code below performs hierarchical clustering analysis on a gene expression dataset and then further analyzes the clusters to identify enriched biological terms using the databases: Gene Ontology (GO), with three Sub-Ontologies (Biological Process (BP), Cellular Component (CC), Molecular Function (MF)) transcription factors (TF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.

Code 
# Get the cluster assignments
mt <- as.hclust(heatmap$rowDendrogram)

# Cut the tree into 8 clusters
tgcluster <- cutree(mt, k = 8)
tgdegnames <- rownames(DEGsFrame)

# Keep unique cluster numbers
cl <- as.numeric(names(table(tgcluster)))

totalresults <- 0
totalcols <- 0
pcols<-c("firebrick4", "red", "dark orange", "gold","dark green", "dodgerblue", "blue", "magenta", "darkorchid4")

# Iterating through clusters for functional enrichment 
#  The gost() function queries the significant genes from each cluster against various biological databases
for (i in c(6, 5, 4, 3, 2, 7, 1)) {
  gobp <- gost(query = as.character(tgdegnames[which(tgcluster == cl[i])]), organism = "mmusculus", significant = T, sources = "GO:BP")$result
  gomf <- gost(query = as.character(tgdegnames[which(tgcluster == cl[i])]), organism = "mmusculus", significant = T, sources = "GO:MF")$result
  gocc <- gost(query = as.character(tgdegnames[which(tgcluster == cl[i])]), organism = "mmusculus", significant = T, sources = "GO:CC")$result
  tf   <- gost(query = as.character(tgdegnames[which(tgcluster == cl[i])]), organism = "mmusculus", significant = T, sources = "TF")$result
  kegg <- gost(query = as.character(tgdegnames[which(tgcluster == cl[i])]), organism = "mmusculus", significant = T, sources = "KEGG")$result
  
  # Combine results
  results <- rbind(kegg, tf, gobp, gomf, gocc)
  
  # Filter the results based on different sources
  tf  <- grep("TF:", results$term_id)
  go  <- grep("GO:", results$term_id)
  kegg<-grep("KEGG:", results$term_id)
  
  # Get enriched terms/pathways, their associated p-values, and other relevant information obtained from the enrichment analysis.
  kegg<- results[kegg, ]
  tf  <- results[tf, ]
  go  <- results[go, ]

  # Order the results based on p-values
  kegg<- kegg[order(kegg$p_value), ]
  go  <- go[order(go$p_value), ]
  tf  <- tf[order(tf$p_value), ]
  
  # Split the term_id and term_name
  ll <- strsplit(as.character(tf$term_name), ": ")
  ll <- sapply(ll, "[[", 2)
  ll <- strsplit(as.character(ll), ";")
  tf$term_name <- sapply(ll, "[[", 1)
  
  # Remove duplicates
  if (length(tf$term_id) > 0) {
    uniqtf <- unique(tf$term_name)
    tfout <- 0
    for (ik in 1:length(uniqtf)) {
      nn <- which(as.character(tf$term_name) == as.character(uniqtf[ik]))
      tfn <- tf[nn, ]
      inn <- which(tfn$p_value == min(tfn$p_value))
      tfout <- rbind(tfout, head(tfn[inn, ], 1))
    }
    tf <- tfout[2:length(tfout[, 1]), ]
  }
  
  results <- rbind(head(kegg, 10), head(go, 10), head(tf, 10))
  totalresults <- rbind(totalresults, results)
  n <- length(results$term_id)
  totalcols <- c(totalcols, rep(pcols[i], n))
}

totalresults <- totalresults[2:length(totalresults[, 1]), ]
totalcols <- totalcols[2:length(totalcols)]
par(mar = c(5, 15, 1, 2))
Code 
# Visualization of under-expressed clusters of DEGs
barplot(
  rev(-log10(totalresults$p_value[75:126])),
  xlab = "-log10(p-value)",
  ylab = "",
  cex.main = 1.3,
  cex.lab = 0.9,
  cex.axis = 0.9,
  main = "Under-Expressed Clusters",
  col = rev(totalcols[75:126]),
  horiz = T,
  names = rev(totalresults$term_name[75:126]),
  las = 1,
  cex.names = 0.6
)
Figure 13: Under-expressed clusters of differentially expressed genes after functional enrichment analysis
Code 
# Custom layout / margins
par(mar = c(5, 15, 1, 2))

# Visualization of over-expressed clusters of DEGs
barplot(
  rev(-log10(totalresults$p_value[1:74])),
  xlab = "-log10(p-value)",
  ylab = "",
  cex.main = 1.3,
  cex.lab = 0.9,
  cex.axis = 0.9,
  main = "Over-Expressed Clusters",
  col = rev(totalcols[1:74]),
  horiz = T,
  names = rev(totalresults$term_name[1:74]),
  las = 1,
  cex.names = 0.6
)

References

Karagianni, Niki, Ksanthi Kranidioti, Nikolaos Fikas, Maria Tsochatzidou, Panagiotis Chouvardas, Maria C. Denis, George Kollias, and Christoforos Nikolaou. 2019. “An Integrative Transcriptome Analysis Framework for Drug Efficacy and Similarity Reveals Drug-Specific Signatures of Anti-TNF Treatment in a Mouse Model of Inflammatory Polyarthritis.” Edited by Cinzia Cantacessi. PLOS Computational Biology 15 (5): e1006933. https://doi.org/10.1371/journal.pcbi.1006933.
Monaco, Claudia, Jagdeep Nanchahal, Peter Taylor, and Marc Feldmann. 2014. “Anti-TNF Therapy: Past, Present and Future.” International Immunology 27 (1): 55–62. https://doi.org/10.1093/intimm/dxu102.
Morgenthaler, Stephan. 2009. “Exploratory Data Analysis.” WIREs Computational Statistics 1 (1): 33–44. https://doi.org/10.1002/wics.2.